""" A module of stuff for doing Shepard interpolation.

Features:
  - Inverse distance interpolation (DONE)
  - Cost cutting by ranking or radius search (DONE)
  - Direction and masking (...)
  - ...

See the Shepard paper.

The function you probably want is L{interpolate()<interpolate>}

  """

from random import random
import sys

#Shepard interpolation

import Image, ImageDraw
from Scientific.Geometry import Vector
from math import cos,sin


def distance(p1,p2):
    """ Returns the distance between two points.
    Points are represented as 2D sequences (tuple, list, etc)
    @type p1: sequence
    @param p1: First point
    @type p2: sequence
    @param p2: Second point
    @rtype: number
    @return: the L2 distance between p1 and p2
    
    """
    return ((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)**0.5


def reduce_points_radius(pos,points,radius): # Cp = {Di|di<=r}
    """Returns a (possibly) shortened list of points from C{points} based on their distance from C{pos}.
    Points that are greater than C{radius} away from C{pos} are discarded.
    If this results in no points, the nearest point is added to the returned sequence.
    @param pos: the centre of the region of interest
    @type pos: sequence (x,y)
    @type points: squence of points ((x,y,v),(x,y,v),...) where v is the value at that point
    @param points: the points to select from
    @type radius: number
    @param radius: maximum distance for inclusion
    @rtype: sequence of points
    @return: points within the radius or the single nearest point
    """
    localpoints=[]

    x,y=pos
    for p in points:
        d=distance(pos,p)
        if d<=radius:
            localpoints.append((d,p))

    if len(localpoints)==0:
        nearest=points[0]
        for p in points[1:]:
            if distance(pos,p)<distance(pos,nearest):
                nearest = p
        localpoints.append((distance((x,y),nearest),nearest))


    localpoints.sort()
    return localpoints

def reduce_points_ranking(pos,points,radius): # CpN = {Di,D1,...,Din} (n<=N)
    """Returns a (possibly) shortened list of points from C{points} based on ranked distance from C{pos}.
    Only the C{radius} number of points closest to C{pos} are returned.  If C{radius} is 0, all points are returned.

    @param pos: the centre of the region of interest
    @type pos: sequence (x,y)
    @type points: squence of points ((x,y, v),(x,y, v),...) where v is the value at the point
    @param points: the points to select from
    @type radius: number
    @param radius: number of points to select
    @rtype: sequence of points
    @return: points selected as nearest
    """

    neighbours=radius
    x,y=pos
    if neighbours==0: neighbours=len(points)
    localpoints=[]
    for p in points:
        localpoints.append((distance((x,y),p),p))
    localpoints.sort()
            
    return localpoints[:neighbours+1] # So that we can get r'(CpR) = min{dij\Dij does not belong to CpN} = din+1

    
def interpolate(pos,points,expo,radius=100,reduce_points=reduce_points_radius):
    """Finds the interpolated value for point C{pos}, based on the known points C{points}.
    The exponent given affects the falloff of effect of each point.  High exponents (e.g. 70) lead to hard edges between regions around a point.  Low exponents (e.g. 2) lead to soft edges.
    
    Optionally, the radius of effect can be given.  If this is set to 0, only ther nearest point is used.
    
    Also, the method of discarding points to affect the interpolation can be changed by changing the default value of C{reduce_points} from L{reduce_points_radius} to L{reduce_points_ranking} or other suitable function.  Note that this also changes the meaning of C{radius}.  See the documentation for L{reduce_points_ranking} for more information.

    @param pos: position of the required interpolated value
    @param points: list of known points
    @param expo: exponent for falloff
    @param radius: radius of inclusion for known points
    @param reduce_points: function to be used to limiting the included points
    @return: interpolated value
    """
#    print "!",radius
    x,y=pos
    val1=0.0
    val2=0.0
    num = 0.0
    den = 0.0
    col=-1
    neighbours = radius
    radius = (7*((256*256)/3.14159*len(points)))**0.5

    sval = [] #Keep calculated si values (di,(xi,yi,zi),si)
    tval = [] #keep calculated ti values ((xi,yi,zi),ti)
    
    cp = reduce_points_radius(pos,points,radius)
    cpn = reduce_points_ranking(pos,points,neighbours)
    
    #cpp = C'p
    #rpp = r'p
    if 0 <= len(cp) <= 4:
        cpp = cpn[:4]
        rpp = cpn[4][0]
    elif 4 < len(cp) <= 10:
        cpp = cpn[:-1]
        rpp = radius
    elif 10 < len(cp):
        cpp = cpn[:10]
        rpp = cpn[10][0]

            
    for p in cpp:

        di=float(p[0])#distance((x,y),p)
        p1=p[1]                
        if di==0: #On a known point
            col=p1[2]
            break
        if 0 < di <= (rpp/3):
            si = 1/di
        elif (rpp/3) < di <= rpp:
            si = (27/(4*rpp))*((di/rpp)-1)**2
        elif rpp < di:
            si = 0
        sval.append((di,p1,si))
    
    for t in sval:
        Di = float(t[0])
        Xi,Yi,Zi = t[1]
        si = t[2]
        for a in sval:
            dj = float(a[0])
            xj,yj,zj = a[1]
            sj = a[2]
            #if not t[1] == a[1]:
            n1 = sj*(1-(((x-Xi)*(x-xj)+(y-Yi)*(y-yj))/(Di*dj)))
            num+= n1  #numerator
            den+=sj   #denominator
    
        ti = num/den
        wi=(1+ti)*si**2
        tval.append((t[1],wi))
        
    for b in tval:
        xb,yb,zb = b[0]
        v1 = zb*float(b[1])
        val1+=v1
        val2+=b[1]
    
    if col==-1: #If not on a known point
        col=val1/val2
    return col



def mappoints_simple(points,im,expo,neighbours=0,reduce_points=reduce_points_radius):
    """
    Given a set of points, an image and the interpolation parameters, this function performs a point-by-point interpolation and draws the map on the given image.

    Mainly used for demoing, at the moment.


    Note
    ====
    
    Prints out percentage completion.  This should probably only happen when a debug flag is set.

    @param points: list of known points
    @param im: Image to draw on - see the U{Python Imaging Library<http://www.pythonware.com/products/pil/>}
    @param neighbours: distance/number of neighbours included in the calculation
    @param reduce_points: function to use when deciding which points to include in the calculation
    """
    map_draw=ImageDraw.Draw(im)
    for x in range(im.size[0]):
        print "%2d%%" %(100*float(x)/im.size[0]),
        sys.stdout.flush()
        for y in range(im.size[1]):
            col=int(interpolate((x,y),points,expo,neighbours,reduce_points))
            map_draw.point((x,y),fill=(col,col,col))

        print "\b\b\b\b\b",
        sys.stdout.flush()
    print
    
def select(probability):
    """Returns true if a randomly generated number is less than the parameter C{probability}.
    @param probability: desired probability of a True result
    @type probability: number, 0..1
    """
    return random()<probability
        


def simulate_local(points,focus,im,drawnodes=True,neighbours=0,reduce_points=reduce_points_radius):
    """ Simulates (poorly, but good enough for visual demo) how a map would look based on information that gets more sparse at greater distances.  Draws on a given image.  Mainly for demo.
    """
    newpoints=[]
    for p in points:
        #Probability of being counted decreases with distance
        prob=distance(p,focus)
        prob/=(im.size[0]**2+im.size[1]**2)**0.5
        prob=1.0-prob
        #Exponential drop
        prob**=3
        if select(prob):
            newpoints.append(p)            


    mappoints_simple(newpoints,map_im,exponent,neighbours=neighbours,reduce_points=reduce_points_radius)
    map_draw=ImageDraw.Draw(map_im)
    #if drawnodes:
      #draw known points
      #for p in newpoints:#points:
          #map_draw.ellipse(((p[0]-1,p[1]-1),(p[0]+1,p[1]+1)),outline=(255,0,0),fill=(int(p[2]),0,0))